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A new cubature formula for weight functions on 
the disc, with error estimates 

O. Kounchev and H. Render 


Abstract 

We introduce a new type of cubature formula for the evaluation of an integral over the disk 
with respect to a weight function. The method is based on an analysis of the Fourier series of 
the weight function and a reduction of the bivariate integral into an infinite sum of univariate 
integrals. Several experimental results show that the accuracy of the method is superior to 
standard cubature formula on the disk. Error estimates provide the theoretical basis for the 
good performance of the new algorithm. 

1 Introduction 

Recently, methods for the numerical evaluation of integrals of the form 

p /' 27 T nR 

R ( 5 ) = / 9 {x) dx = / g {r cos (p,r simp) rdrd(p (1) 

Jd Jo Jo 

on the disc Dn of radius R in the plane have received increased attention in the framework of the 
meshless local Petrov-Galerkin (MLPG) method, see [15], [T^], [3S], [5S] . Numerical experiments 
in m have given evidence that classical rules like the piecewise midpoint quadrature rule, or the 
rule of Peirce (for definitions see below (l3^ and (l40l) l are superior to the Gauss-Legendre product 
rule which is very popular in the MLPG literature. 

In the present paper we study new methods for the numerical evaluation of integrals of the 
type 

Iw if)=[ f (x) w (x) dx (2) 

JDr 

where w {x) is a (not necessarily non-negative) weight function on the disc Dr in the plane The 
introduction of a weight function is an important concept in numerical integration: the integrand 
g (x) is decomposed into a product / (x) w {x) where the factor / (x) is a function well-approximable 
by polynomials (see section 3.7 in [46) 1 and w (x) is a function of limited smoothness or with a 
singularity. Using the specific properties of the weight function w (x) one aims to achieve a cubature 
formula for the integration of the function / with respect to w (x) dx which should be more accurate 
than using directly a cubature formulae for g = f ■ w like in O- 

The main concept underlying our method is to expand the weight function w (x) = w (re*"^) into 
a Fourier series o and rely upon a similar expansion for the polynomial-like function / (x), called 
the Almansi expansion, see ®. The reader will find an explicit description of this construction 
below, after all necessary notations and tools are introduced. Illustrating examples in this paper 
are the weight functions 

{x,y) = and {x,y) = \y\ (3) 


1 




where the first weight function has a singularity in 0 and the second weight function is continuous 
on the closed ball but is not differentiable on the line ?/ = 0 in the interior of the disk. We shall show 
by numerical experiments, and by theoretical considerations as well, that our method provides in 
these cases results which are better than the above-mentioned methods for the approximation of 
the integral for the function g (x) = f (x) w (x). 

Let us now give a detailed introduction to the main topic of the present paper, the numerical 
evaluation of the integral with respect to a weight function w. A central role in our approach plays 
the Fourier series of w, given by 

oo as, 

w (re*"^) =w(r cos (p, r sin ip) := (r) {<p) ■ (4) 

fc=0^=l 

Here we use a notation for the Fourier series which is more convenient in our context, and which is 
well known from the theory of spherical harmonics (these convenient notations will be important 
also for further multivariate generalizations, as in [30)1: we work with the orthonormalization of 
the harmonics cos kp and sin kp, defined by 

y(o.i) {p) = i/V^ (5) 

Y(ku {p) =-^coskp and Y/^p) (p) =-^ sin kp. (6) 

VtT y'TT 

for integers k > 1. Then Y(^kp-) is an orthonormal system for fc > 0,£ = l,..,afe, where = 2 
for fc > 1, and oq = 1. The {k,i)-th Fourier coefficient of a complex-valued continuous function 
/ (re**^) is 

p27r 

fik,i){r):= f {re’‘‘^) Y(^k/) {‘P) dp for fc > 0, £ = 1,Ofc (7) 

^0 

and the corresponding Fourier series of / is 

oo flfc 

/ = fir cos p, r sin p) := iO (p) ■ (8) 

fe=0^=l 

Let us recall that the Fourier series of a polynomial p (x, y) is of a very special form: there exist 
polynomials P(k,e) and a number N < degp{x,y) such that 

N aji, 

p (x, y) =pir cos p, r sin p) = ip ); (9) 

k=0t=l 

the representation (jH]) is called Gauss decomposition or Almansi expansion of a polynomial p. 
Hence, the Fourier coefficient P{k,t) (r) of a polynomial p (x, y) is of the form 

P{k,l) (r) = P{kp) (r^) r'". (10) 

Moreover the degrees of all P(k,e) are bounded by — 1 if and only if the polynomial p (x, y) is 
polyharmonic of order TV, i.e. if A^p {x, y) = 0, where is the N-th. iterate of the Laplace 
operator A. These results are given a thorough treatment in see [27], m, m- 

We consider now the integral ([2]), which after introducing polar coordinates, becomes 


p2T7 pR 

Iw{f)= / / f [rcosp^rsinp) ■ w{rcosp,rsinp) ■ rdrdp. 

Jo JQ 

We replace the weight function w (re*‘^) by its Fourier series, and after interchanging summation 
and integration we obtain 


/»(/) 



oo gfc 

w (re*"^) rdrdp = FF / iO '^{k/) (r) rdr. 

k=oe=i 


( 11 ) 
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Note that for the constant weight function w {x^y) = 1 = one obtains simply 

/•27r pR pR 

hif)= / / f {r cos ip, r sin ip) r dr dip = /(o,i) (12) 

Jo Jo Jo 

Formula (EH) is central to our approach since it reduces the integration of the bivariate function 
/ (re*"^) with respect to a weight function w to the calculation of an infinite family of 

univariate integrals with weight functions W(^k,i) {’>')■ 

Here we come to the most crucial point of our approach: as we assume that the function / 
is well-approximable by a polynomial p, we know that 1^ (/) is close to lyj (p). By the Almansi 
formula m and EHD, the one-dimensional integrals in (fTTl) for computing {p) are equal to 


pR pR 

/ Pk,e {r) W(^k,i) {r) rdr = / pk,e {r'^) r^wt^u,e.) {r) rdr. 

Jo Jo 

After a change of the variable p = ii becomes 

r2 

^ Jo P{P)- p''^‘^W(k,l)[y/p)dp. (13) 

Now we want to employ the iV-point Gauss-Jacobi quadrature to the last integral, with measure 
(\/p) ■ 1'®'^ reason we need the following assumption which will be made throughout 
the entire paper and which is called the pseudo-definiteness of the weight function: 


General Assumption: Each Fourier coefficient W(^k,e) (f) of the weight function w does not 
change the sign over the interval (0, R), and it is integrable and continuous over (0, R). 


Due to our assumption we infer the existence of the A-point Gauss-Jacobi quadrature with 
nodes and coefficients (which are either all positive or all negative) 

^l,(fc.£) < (14) 

(15) 

Due to the exactness of the Gauss-Jacobi quadrature for any integer 0 < s < 2A — 1 we obtain 
the equalities 


• ^,{k,i) = 9 / p''p''^^w^k,e) iffip) dp= r^^r'^unk,i) (r) rdr. (16) 

^ Jo Jo 

Hence, for a polynomial /, for which deg/(fe^^) < 2A — 1 for all {k,i), by using the Gauss-Jacobi 
quadrature EH), the integral EH becomes 

Iw if) = / hkj) (^) (?') rdr (17) 

k=oe=i 

, oo ak 

= 2 ^^ / iVp)}dp 

k=0i=l 

= U) 

where we have put 

In'"' (/) := ' h,{k,e) ' f(kJ) iV^Jk/)) ■ (18) 

k=oe=ij=i 
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In [20], we defined 1^°'^ (/) as polyharmonic cubature of degree N in arbitrary space dimension. 
The reason for the name is the fact that is exact on the space of all polynomials of polyhar¬ 
monic order < 2N, i.e. for each polynomial / such that = 0, or more general, on the space 

of smooth functions satisfying the polyharmonic equation = 0 

In previous work [28], [29] . m, we have given a motivation and a detailed analysis of the 

polyharmonic cubature of degree 2N in the framework of the Polyharmonic Paradigm, explaining 

_^ 

the natural appearance of the factor t^ in our formulae which is related to the Gauss-Almansi 
decomposition of a polynomial. 

As the values f{k,£) s-re Fourier coefficients, they can be approximated by means of 

Discrete Cosine/Sine transform of the function /, by which we mean the expression 


AM) 


(r) 


n M 


Y(k,i) 



(19) 


The main contribution of the present paper is the Discrete Polyharmonic Cubature with 
parameters {N, M, K ), defined for integers iV > 1, M > 1, and Ai > 0, by putting 


K ak N 


(/) := • h,h) ■ /(M) (vE^) 

K ak N M 


TT 

M 


a k ivi / O \ 

EEEE^J.(M) • *j,(k,e) ■ Y(k,e) {-Jfj • / ■ 


k— 0 £—lj—ls—l 


( 20 ) 

( 21 ) 


We see that unlike (fTSl) formula (1^ for (/) is indeed a cubature formula in 

sense of the word. Its coefficients (weights) • t. • Y(^k,e) (tt) } have varying 

they satisfy the following remarkable inequality 


TT 

M 


K ak N M 

EEEE 

k—O^—lj—ls—l 




_ k 

■ ■ Y{k,e) 



< ^Ikll, 


the usual 
signs but 


( 22 ) 


where it is assumed that the weight function w satisfies the so-called summability condition 

OO ak pR 

11 ^11'= EE/ \w(k,e)ir)\rdr <oo. (23) 

k=0 £=1 

Hence, if the weight w satisfies the summability condition, the cubature coefficients satisfy the 
important stability inequality (1221) . and experimental evidences show that all coefficients are in 
general very small. This inequality is proved in Theorem [2| by an application of the famous 
Chebyshev extremal property for the Gauss-Jacobi quadrature, see Theorem 4.1 in Chapter 4 of 

[33]. 

Let us give a short outline of the paper. 

In Section [2] we provide basic properties of the discrete polyharmonic cubature formulas: the 
summability condition for the weight function implies that are uniformly bounded func¬ 

tionals (in the parameters N, M, K) on the set of all polynomials. Under this assumption it follows 
that for N,M,K —oo, the value k) (/) converges to (/) for any function / which is 

continuous on the closed disk with radius R. Moreover we show in Theorem |4| that our cubature for¬ 
mula X) exact for all polynomials of the type / (x) = r'^‘^^^Y(^k,e) (,a) where 0 < s < 2N— 1, 
0<k<'M—1 — K, and £ = 1,2, ..., ak, i.e. 
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Section [3] is devoted to error estimates for the discrete polyharmonic cubature. The error 
bounds are a sum of the error bounds of three successive approximations: 1. the approximation 
of the weight function w as a Fourier series - involving the parameter K; 2. the approximation by 
the one-dimensional quadrature formula in radial direction - involving the parameter N; 3. the 
approximation by the Discrete Fourier Transform - involving the parameter M. 

In Section|3]we provide experimental results for the discrete polyharmonic cubature with respect 
to the first weight function in (|3]) for four different types of test functions, and compare the 
results with those obtained by the piece-wise midpoint rule ([H]) and the rule of Peirce m), 
which are two methods used widely in practice, in particular in the Meshless Petrov-Galerkin 
method, see m, ca, m, isa- Our methods have much higher accuracy than all other methods 
which might be explained by the fact that the weight function has a discontinuity in 0 which 
has a strong negative influence on the results of the usual cubature formulas. 

Section [5] contains experimental results for the second weight function in ([3]). Since the 
weight function has an infinite Fourier series it has been expected that the discrete polyhar¬ 
monic cubature now would perform weaker than for the first weight function. However, even in 
this case the accuracy has been extremely good provided that the parameter K is large enough. 
Finally, in Section [5] we comment on practical aspects for the numerical implementation. 

2 The discrete polyharmonic cubature 

At first we introduce the notion of pseudo-positive (or more general, pseudo-definite) weight func¬ 
tions which plays a central role in the discrete polyharmonic cubature formula. 


Definition 1 Let f (x, y) be a function on the disc with center 0 and radius R which is continuous 
except for 0. Then f is called pseudo-definite if every Fourier coefficient f(^k,e) (r) of f has a definite 
sign, thus either f(k,e) (r) > 0 for all r G (0, R) or f(^k,i) (r) < 0 for all r G (0, R ). The function f 
is called pseudo-positive if f(k,e) (r) > 0 for all r G (0,i?) for all k G No, I? = 1, ...,afc. 

A simple example of a pseudo-positive function is the Poisson kernel on the unit disk given by 


P{x,y) 


1 — x'^ — y'^ 
{x- 1)^ -f 


1 _ J.2 

1 — 2r cos ip -\- r"^ 


OO 

1 + ^ 2r^ cos kip. 


Let us recall that a cubature rule is just an expression of the type 

N 

a if) = '^cjf{xj) 
i=i 


where xi,...,xn are pairwise different points in K”, called nodes, and Cj ^ 0 are real constants 
called weights. An important property of a cubature rule, which is also a basis for widely used 
method for construction of cubature formulas, is the exactness for a subspace of functions: a 
cubature formula C is exact on a subspace of functions U with respect to a weight function w [x) 
if 


C {f) = J f (x) w (x) dx for all f G U. 


If U is the subspace Vm of all polynomials of degree < m then one says that a cubature formula 
C has degree to if C is exact on Vm and if there exists a polynomial / of degree to -I- 1 such that 
C if) Iw if) ■ Exactness of degree to is also useful for providing error estimates using Taylor 
series, an approach which was emphasized already by R. von Mises, see [38], [39]. 

Recall that the Discrete Polyharmonic Cubature formula with parameters {N,M,K) is 
defined for any continuous function / on the disk of radius R by formula m- 


rpoly 


if) 


^ K ak N M ^ 

k—Ol—lj—ls—1 





(24) 
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Thus if) is a cubature formula where the weights are real numbers. At first we establish 

the following subtle estimate (analogous to an estimate proved for the polyharmonic cubature 
formula in (SE): 

Theorem 2 Let w he a pseudo-definite weight function with Fourier coefficients W(^k,i )7 o,nd let f 
be bounded on the closed disk of radius R with supremum norm 

ll/lloo := sup \fix,y)\. 


Then 


^ pR 


pJX 

^(N%,K)(f) \wik,i)ir)\rdr < ^/f^\\f\\^\\w\\ 

1 _n f) _1 0 


Also, the coefficients of formula ^2^ satisfy inequality 

Proof. Since |T(fc,£) (v?)! < for all tp € [0,27r] the following estimate follows directly from 

K ak N ^ 

ll/lloo I 

fc=0^=l:; = l 
N 

Following (ITSl) we define the function (h) := \ ^j,(k,e) \ h {tj,{k,i)) i which due to the pseudo- 

i=i 

definiteness of w, is the A^-point Gauss-Jacobi quadrature for the integral 


1 

i{k,i) Qi) -= 2 ]^ ^ (\/^)| 

Now we apply the Chebyshev extremal property of the Gaufi-Jacobi quadrature (see Theorem 4.1 
in Chapter 4 of [33]) which shows that 

ih) < I^k/) ih) 

holds for any 2A^—smooth function h (r) with (p) > 0 for all p > 0. Applying this to the 

function h (p) = p~^l'^, gives the inequality 

(p"G2j ^ |Ay(fc,£)| < - / \w(^k,i) iVp)\dp= \w{k,e) (r)| rdr. (25) 

j—1 *^o 


From Theorem[2]by using standard results from functional analysis (see [131 P- 351]) we obtain 
the following important: 

Corollary 3 Suppose that the Fourier coefficients of the pseudo-definite weight function w satisfy 
the summability condition (see \23\) ) 

||rc|| < 00 . 

Then the discrete polyharmonic cubature (/) converges to the integral (/) for each 

continuous function on the closed disc with radius R when N, M, K are approaching infinity. 

^The notations in this reference differ from the present: Here we have put drpk,!,) (r) = f^^ (v) dp (re*'*’) 

while in m we used to work with (r) = ^ik,£) . 
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Let us recall that the trapezoidal sum (see [T71 p. HI]), also called the composite trapezoidal 
rule in [231 p. 155], is defined on the set of all 27r-periodic functions g on the line, by 

, , 27r / 27rs\ ^ ^ 

S = 1 ^ ^ 

Let now f{x,y) be dehned on the disk with radius R and dehne fr (ip) := f (r cos (p,r sin ip) for 
r G [0, i?]. Then 


Tm {fr ■ Y^k,i)) = (^cos ,rsin 


2tts 

W 


(^) 


and the discrete polyharmonic cubature (1201) can be written as 

K ak N 


1 k / \ 


fc=o^=ii=i 


For the following result recall that r^Y(^k P) {ip) is a polynomial of degree k in the variables 
X = r cos p and y = r sin p, cf. [27], m- 


Theorem 4 Let M and K be natural numbers satisfying M > K. Then the discrete polyharmonic 
cubature with parameters {N, M, K) is exact on the linear subspace generated by the 

polynomials 

r^^+'^Y^k,e) (t) 

where Q < s <2N — 1, k < M — 1 — K, and £ = 1, ...,ak- 


Proof. We shall use that Tm {g) = g (p) dp for all trigonometric polynomials g of degree 
k < M — 1, see e.g. HZ], p. 110. Let k < K he a. natural number and let fci < M — 1 — K he 
a natural number. Then it is easy to see that the product (p) (p) is a trigonometric 

polynomial of degree fci + fc < M — 1. It follows that 

p2-K 

^M,0 (^(A:i/i)^(fc/)) = / (28) 

^0 

where Stj is the Kronecker symbol. Now we consider the polynomial 

fisM,ei) ^rcosp,rsmp) = r^^r’^^Y^^.A) {t) 
with s < 2N — 1 and ki < M — 1 — K. Then by (1^51) we obtain 


rpoiy 

\n,M,K) 


K ak N 


^ {tjpk,e) {\/^ 3 ,{k/)) Y{kiA)'^ik,i)) 

k=oe=ij=i 

^ ^ p27r 

= / ^(fci.G) (t) Y(kA (t) dp 

^ k=oe=ij=i 

1 ^ 

j=i 

= ^ (^) "^dr = flkf}f)"^ (r) (r) rdr 
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In the last row we used formula (HU), due to the fact that the Gauss-Jacobi quadrature for each 
(fci,£i) is exact for all polynomials of degree 2N — 1. In the Fourier series expansion of the polyno¬ 
mial the only non-zero coefficient is {r) = Hence, by formulas (ITT)) . ((T5)) 

we obtain 

rpoly f f(s,ki,£i)'\ _ j ( 

^{N,M,K) J — yj J ■ 


Remark 5 In one can find a nice account of the early history of numerical multivariate inte¬ 
gration until the 1950’s, commencing with the work of J.C. Maxwell f37| / in 1877, of P. Appel B 
m and H. Bourget m- A. Ahlin m stresses the fact that for many ad hoc formulae there are no 
error estimates available, and he provided error estimates for the case of product type measures. 
Since the 1950’s the literature has grown considerably and very good surveys until the 1970’s can 
be found in the books UB’ ® UB and m- For a recent survey on numerical integration 
rules we refer to m- For numerical integration rules especially for the disc we refer to UB, UM, 
m, m„ m and older work in m, m , Wi see also JTB> M- 

3 Error Estimates for the discrete polyharmonic cubature 

In this section we derive error estimates for the discrete polyharmonic cubature by considering 
the following chain of identities and approximations: at first the integral ([U is equal to (fTTl) . and 
a simple change of variables leads to (1291) : then the infinite sum is approximated by a finite sum 
(depending on K), and next we employ the Gauss-Jacobi quadratures (depending on N), and 
finally we use the ’’Discrete Fourier transform” (depending on M) in order to obtain the discrete 
polyharmonic cubature (/): 


'Dr 


f (x) w (x) dx 


PTt ^ nJX 

= XI / AfcJ) (^) (^) = o X X / 

k,i ^ k^O 

^ K (Ik 

f(k,e) (Vp) P " X p'^w{k,e) (Vp) ^P 

^ k=0 £=\ 

if afc AT _ ^ 

^ 2 ^ ^ ^ ^ (vX(M)") ^ ^{N^oo,K) (/) 

k=0 i=l j=l 

.. K ak N ^ 

~ 2 X X X Am) ^ = ^(nJm,k) (/) 


k=oe=ij=i 

The error between formulas (ED and (1321) can be estimated as follows: 
Theorem 6 For f G 0 "^^+^ {Br) the following estimate holds 

2 ttC{ 2 D + 1 )^ 


jpoiy ( f \ — rP°'y ( 
^(N.oo.K)\J i > 


< 


M'^d +1 


XX^'^-^ (/) / W(k,t) {r) rdr 

1 -n D _1 0 




where 


Mk,i if) ■■= sup sup 

re[0,i?] (/3€[0,27r] 


^ 2 D+l 


dif'^D+l 


[f (re*"^) Y(fc,£) (v?)] 


Here f denotes the Riemann zeta function. 


(29) 

(30) 

(31) 

(32) 
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Proof. In [T71 p. 110], for periodic g G ( 720+1 27 r] it is shown that 


^271- 


9 iv) d(p - Tm (g) 


< 


a 


g,D 


jVf 20+1 ’ 


(33) 


where the constant Cg^D is given by 


Cg^D := 47rC {2D + 1) sup 

c/?G[0,27r] 


^ 20+1 


[9 (</^)] 


Formula (|27|l and the triangle inequality show that for = / (re®'^) 


K ak N 


'^(Ar,L,K) (/) (/) — (\/^i.(fe7)) 

fc=0^=lj=l 

Now we use inequality (1331) for estimating the term in the middle row: 


< 


4< (2D + 1) (/) 


^ 20+1 

Hence, by inequality (1^ . we obtain the final result: 


rpoly 


rpoiy 


, (f) — P 

‘(Af.oo.iC) w ^{N,M,K) 


if) 


27rC {2D + 1) j 


< 


M-^d+ 

2'kC {2D + 1) 

M^o+i 


fc=0^=l 
K ak 


i=i 

rR 


^k pn 

(/) / lU(fcy) (r) rdr. 
k=oe=i do 


The error bound between (|29l) and (l30)l is considered in the following: 

Theorem 7 Let f G (D/j) for some integer p > 1. Assume that the weight function w satisfies 

^ 00 afc ./{S 00 Ofc 

/ h(M) (v^) I = XI X! / h(fc 7 ) = ||u>|| < 00 

7_n D _n do 7-n D _n dO 


fe=0 ^=1 


fe=0 ^=1 ■ 


Then 


00 Qfc 

E E /(fe7) (Tp) P = X p'^W(^k,e) (Vp) dp 
, 7 «_7 do 


< 


1 ^ 






fc=K+l ^=1' 

Proof, By applying a standard techniques (see e.g. Theorem 10.19 in m), we obtain 

|.27r 1 |.27I- ^2p 

/(fc7) {r) = f (re**^) 4"(fey) (<p) dp=^ f (re*"^) ^Yik.i) (^) d<P 


1 

k2p 


d^P 

dp^p 


f {re^'P) Y(^k,e) {p) dp- 


We can estimate this inequality by using the Cauchy-Schwarz inequality and the orthonormality 
of Y(k,e.) (<p), arriving at 


If ( W ^ 
|/(fe,f) (r)| < 


E. 

dg? 


f (re^^) 
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Then 


OO gfc pR 

E E/ f{k,t) iyfp) Wik,e) (Vp) M 

=K+1 1=1 
OO ak 

^ E E/ \fik,i) {Vp) \ \w{k,t) {\fp)\ dp 

1 _ \ 1 n _1 J 0 


k ^ K+1 ^=1 


< 


< 


dip' 


& 




OO ak 1 

S J2w, 


rR 


< 

This ends the proof. 




^ ^ \'^{k,i){^fp)\dp 

°° k=K+li=l 

°0 °-k 


1 


X 


if2p 


E E/ \w(k,e) {Vp)\ dp 

— /’—I 






The error bound between (1501) and (I3ip is a Markov type estimate (see Theorem 14.2.2 in m, 
and Theorem 44 in ED). Using the notations (IT4l) . (IT5|) . (ITbl) we prove: 

Theorem 8 Let KN,{k,i) be the leading coefficient of the N-th degree orthonormal polynomial 
QN,(k,t) with respect to the measure p^W(^k,e) (\/p) Ibe interval [0,i?^]. Then for every function 
f G {dJR) and every index {k, £) we have 


‘-(k/) 


N 


< 


f{k,i) i^fp) W(k,i) (Vp) dp-"^ f{k,i) {i/tj,ik,e)) X ijXk,e)^j,ik,i) 

i=i 

d^^ 


1 


(2iV)! 


Ik?. 


''N,(k,e) 


dp 


,2N 


9(k,e.) (p) 


where g[k/) (p) '■= p ^ f(k,e) {\fp) ■ The difference between L30\) and L31\) is bounded by 

d^^ 


K ak .. if Ofc 

E E hk,e) < E E 


dp- 


2N 


9{k,e) (p) 


k=o e=i '• ~'-''k=oe=i N,(k,e) 

Proof. It is easy to see that g(k,e) (p) = P~^ f(k,e.) {\fp) is 2iV-times differentiable in the interval 
(O, R^) . Hence, to the function g[k,i) (p) = f{k,e) (\/p) P~^ £ (O; di^) ) may apply Markov’s 

Theorem 14.2.2 in m, and we obtain 

d^^ 


1 

- i2Ny.K% 


dp 


2N 


9{k,i) (p) 


^N,{k,t) 

This implies the estimate for the error between (1501) and (15T1) 

K ak .R^ K ak N 


’-^k pJri 2v 

EE/ f(k,l) (p) W(fc.£) (Vp) dp - EEE f{k,i) {i/tj,(k,e)) X tj^(k,e)^3,{k,e) 

I. 1 . n /} 1 _ 1 


k=o e=i 


K ak 

^EE 

k =0 i=\ 


fc=0 1=1 j=l 
N 


fR^ ^ 

9{k,i) (p) X W(^k,i) (f/p) dp ^ ) f{k,t) X tj^(^k,£)^j3k-/) 


i=i 


K ak 


< 


(21V)! 


EE;^ 

k =0 £=1 ^N,{k,l) 


d^^ 

jN9(k,£) (p) 


dp- 





















































Remark 9 In the case when (Vp) Jacobi weight functions explicit expressions for 

KN,{k,i) o^re well known, see e.g. JJ5] / or Theorem 44 in 131}/ . 

Remark 10 If g(k,i) i^) on the real axis and 2 tt periodic and if g(k,e) + w) holomorphic 

in a strip \y\ < a, then one may use estimates in p. 110], and apply them to obtain an 
alternative estimate in Theorem\^ 


4 Experimental results for the weight function (x, y) 

In this section we want to test cubature formulas for integrals of the type 

p2-K nR 

Iw{f)= / / f {r cos g},r simp) ■ {r cos (p,r simp) ■ r dr dip 

Jo Jo 

for the weight function 


{x,y) = 


1 + X 

\/ a ;2 + 2/2 


= - + cost/j = ^^^y(o,i) {p) + V^Yiip) (t) ■ 


Since has only two Fourier coefficients we take iF = 1 in the discrete polyharmonic cubature 
(l20l) with parameters {N,M,K). For given N one has to construct the iV-point Gauss-Jacobi 
quadratures for the integrals 

1 / P {p) p~^^^dp for fc = 0, and ^ / P {p) P^^^dp for /c = 1 

2 7o 2 Jo 

which are Jacobi-weight functions. Fortunately, there are excellent programmes for determining 
the nodes and weights of the Gauss-Jacobi quadrature providing high accuracy, see 
For our experiments in this section we consider four test functions: 


fo ix,y) = l + x‘^ + y^, 
fi (x, y) = 1 + 


2;2 -t y‘ 

f 2 {x,y) = cos (lOx-F 20?/), 
h (x, y) = (x^ + y^f^^ = 


= 1 + cos^ ip -\-r^ sin' ip, 


(34) 

(35) 

(36) 

(37) 


The first test function /o is a polynomial of degree 4, and the second /i is not smooth at 0. The 
function /2 is of oscillatory type and an example of a test function used in the package of Genz, 
see P3], [33] of the form (in our case u = 0) 


cos (27rM -I- ax -I- by) = cos (27ru) cos (ax + by) — sin (27rit) sin (ax -I- by). 


At first we present the experiments for the Discrete Polyharmonic Gubature, and then we compare 
it with two other standard rules which are applied to the functions 

gj (re*‘^) := fj (re*”^) (re*"^) 

Note that gj is not continuous at 0 for j = 0,1,2, and one might argue that the standard rules 
do not perform too well for functions with a discontinuity. For this reason we have included the 
function /a {x,y) for which gs (r'“^) = r^/2 -|- rcosp) is clearly continuously differentiable, and 

our discrete polyharmonic cubature formula m performs better than the usual methods as well. 

Let us note that the discrete polyharmonic cubature formula (1211) needs (at most) (2K — 1) • 
N ■ M evaluation points. Since has only two non-zero Fourier coefficients we need in this case 
at most 2N ■ M evaluations points. 








4.1 Results for the Discrete Polyharmonic Cubature Formula 

In the following tables we present experimental results for the discrete polyharmonic cubature 
where the number M of points on the circles is chosen to be equal to 9, 25,63, 83. The reason for 
this choice is that we used the Fast Fourier transform for the implementation of the trapezoidal rule 
(depending on M). The number N of concentric circles is chosen to be equal to 10,15, 25,35, 50. 

Due to the exactness of the discrete polyharmonic cubature, the value (/o)j according 

to Theorem |4l must be identical with the true value of the integral if M satisfies M > 6, and 
N >2. This is numerically confirmed by our experiments: for all N and M as above specified we 
obtained up to double precision that 

h ~ 6.754424205218060. 

The second test function can also be integrated explicitly and the true value is 

h = —TTRi 6.872 233 92972767. 

Our experimental results are contained in the following table: 


Table 1 


N/M 

9 

25 

63 

83 

10 

6.87224296287783 

6.87224296287783 

6.87224296287783 

6.87224296287783 

15 

6.87223588060173 

6.87223588060173 

6.87223588060173 

6.87223588060173 

25 

6.87223420205342 

6.87223420205342 

6.87223420205342 

6.87223420205342 

35 

6.87223400297000 

6.87223400297000 

6.87223400297000 

6.87223400297000 

50 

6.87223394775545 

6.87223394775545 

6.87223394775545 

6.87223394775545 

Error 





N/M 

9 

25 

63 

83 

10 

0.00000903315016 

0.00000903315016 

0.00000903315016 

0.00000903315016 

15 

0.00000195087406 

0.00000195087406 

0.00000195087406 

0.00000195087406 

25 

0.00000027232575 

0.00000027232575 

0.00000027232575 

0.00000027232575 

35 

0.00000007324233 

0.00000007324233 

0.00000007324233 

0.00000007324233 

50 

0.00000001802778 

0.00000001802778 

0.00000001802778 

0.00000001802778 


Discrete polyharmonic cubature for = (l + cos^ (p + sin^ ip) ^—h 


cos p 


Obviously we obtain very good approximations: even for the case of 180 = 2 x 9 x 10 evaluations 
points the error is smaller than 10“®. Note further that in the rows we do not obtain improvements 
when M is getting larger. This is due to the fact that p i —> fi is a trigonometric polynomial 

of degree < 7, and the exactness of the trapezoidal rule (E51) implies that there is no change when 
M is larger than 7. 

Now we want to test the function ((551) 


h {x, y) = cos {ax + by). 


It is not difficult to see that the first Fourier coefficient of /2 satisfies 

/ 2 ,(o,i) {r) = V^Jo + b^r'j , and 

/2,(1,1) (j’) = /2,(1.2) (j’) = 0, 

where J„ is the Bessel function of the first kind of order n defined by 


Jr 




CXD 

E 

s=0 


(- 1 )^ 


s! (n + s)! 



(38) 


12 





















cf. [S]. Then the first Fourier coefficient 5(o,i) (r) of the function g ■.= f 2 ■ can be computed 

(using the orthogonality relations of spherical harmonics) 


ff(o,i) W = W^«(dli) = Jo (\/a2 + &v) 

It follows that 

Ii = V^TT j 3 ( 0 , 1 ) {r)rdr = J Jo + b'^r^ dr. 

Using the power series expansion of the Bessel function in (1551) the integral can be evaluated up to 
any accuracy. For a = 10 and b = 20 we see that up to double precision we have 

4(1) (/2) = /i = 0.301310 995335 215 

Our experiments for the discrete polyharmonic cubature provide the following table: 


Table 2 


N\M 

9 

25 

63 

83 

10 

-0.08102057453745 

0.31409913156633 

0.30131093100867 

0.30131093100867 

15 

-0.08102397430499 

0.31409919589293 

0.30131099533522 

0.30131099533522 

25 

-0.08102401217317 

0.31409919589293 

0.30131099533522 

0.30131099533522 

35 

-0.08102401237119 

0.31409919589293 

0.30131099533522 

0.30131099533522 

50 

-0.08102401237809 

0.31409919589293 

0.30131099533522 

0.30131099533522 

Error 





N\M 

9 

25 

63 

83 

10 

0.38233156987266 

0.01278813623111 

0.00000006432655 

0.00000006432655 

15 

0.38233496964021 

0.01278820055772 

0.00000000000000 

0.00000000000000 

25 

0.38233500750838 

0.01278820055772 

0.00000000000000 

0.00000000000000 

35 

0.38233500770641 

0.01278820055772 

0.00000000000000 

0.00000000000000 

50 

0.38233500771330 

0.01278820055772 

0.00000000000000 

0.00000000000000 


Discrete polyharmonic cubature for f 2 wJ'^ = (cos (lOx + 20?/)) ( —|-cos((9 


Note that the formula is sensitive with respect to the values of M : if M is 9 then large 
deviations occur (even if N is large) , for M = 25 the approximation error is 0.01 and the number 
N of circles does not influence much the results. For M = 63 and iV = 10 the approximation error 
is very small: 


0.301310 995 335 215 - 0.30131093100867 = 0.000 000 064326 545 


Next we consider the test function /a (x, y) = for which the integrand (1 + r cos tp) 

is smooth. Then the explicit computation gives 




(1 + r cos <p) rdrdip = 1.795 195 802 051 31 


13 























with the following table: 


Table 3 


N/M 

9 

25 

63 

83 

10 

1.79513323182095 

1.79513323182095 

1.79513323182095 

1.79513323182095 

15 

1.79518029482336 

1.79518029482336 

1.79518029482336 

1.79518029482336 

25 

1.79519315318245 

1.79519315318245 

1.79519315318245 

1.79519315318245 

35 

1.79519497859942 

1.79519497859942 

1.79519497859942 

1.79519497859942 

50 

1.79519556405565 

1.79519556405565 

1.79519556405565 

1.79519556405565 

Error 





N/M 

9 

25 

63 

83 

10 

0.000062570230356 

0.000062570230356 

0.000062570230356 

0.000062570230356 

15 

0.000015507227945 

0.000015507227945 

0.000015507227945 

0.000015507227945 

25 

0.000002648868861 

0.000002648868861 

0.000002648868861 

0.000002648868861 

35 

0.000000823451885 

0.000000823451885 

0.000000823451885 

0.000000823451885 

50 

0.000000237995663 

0.000000237995663 

0.000000237995663 

0.000000237995663 


Discrete polyharmonic cubature for ( “ + cosy: 


4.2 Comparison with the piece-wise midpoint rule 

The piecewise midpoint quadrature rule (see e..g. m) is rather geometric: subdivide the disk of 
radius R by concentric circles of radius 


f -j + 1/3 R j 


J - 2 


N 


—R for j = 1, ...,N 


and radial half-lines with angle 2'Ki/M for i = 1, M. Then the integral over each subdomain is 
approximated by the evaluation of the integrand at the centroid of the sector multiplied by the 
area of the sector, given by (j — (^) Formally, we define: 

Definition 11 The piecewise midpoint quadrature rule, is given by: 


rmid 

^N,M 


(/) := 


2ttR^ 
M ■N'^‘ 


N M 


EEU-b / 


• 27r(s-i) 


M 


M 


3 = 1 s=l 

The results for the first test function /o are contained in the following table. 


(39) 


Table 4 


N=M 

Inm 

Error 

5 

6.293948149 525 97 

0.460476055 692 09 

10 

6.552 664 285 742 99 

0.201759919475 07 

20 

6.652 725 619 004 72 

0.101698 586 213 34 

100 

6.733953574 71790 

0.020470630 50016 

200 

6.744180 708 69065 

0.01024349652741 


Midpoint cubature for = (l + ^—h cosy? 

43 

True Value is: — tt r; 6.754 424 205 218 060 
20 


Note that the case N = 200, i.e. 40 000 evaluation points, still leads to an error of 0.01. The 
discrete polyharmonic cubature was exact for this case. 

For the second test function /i we have a similar pattern: 































Table 5 


N=M 


Error 

5 

6.479185 720 36913 

0.393048209358541 

10 

6.671455 800 85980 

0.200778128867871 

20 

6.770 780 856 01381 

0.101453073713858 

100 

6.85177311709146 

0.020460812636194 

200 

6.861992 88760082 

0.010241042126847 


Midpoint cubature for = (l + cos^ </3 + r® sin^ ip) ^—h cos 


True value is « 6.872 233 929 727 67 


and we see that the error is quite big; for 200 x 200 evaluation points it is bigger than 0.01. 
For the third test function /2 (cc, y) = cos (lOx + 20?/) we obtain the following table: 


Table 6 


N\M 9 25 63 83 

10 -0.190440454101284 0.0936727156130806 0.105393884431863 0.105393884431863 

15 -0.120671303989885 0.154856209500921 0.167165382069484 0.167165382069483 

25 -0.0641280979279670 0.207320729809817 0.219935669864855 0.219935669864855 

35 -0.0400423613447862 0.230317338211144 0.243017030992702 0.243017030992701 

50 ^0220967384451548 1^7690217929()60~ 0.260435021528943 0.260435021528943 ~ 

Error 

N\M 9 25 63 83 

10 0.49175144943649 0.207638279722134 0.195917110903352 0.195917110903352 

15 0.421982299325100 0.146454785834294 0.134145613265731 0.134145613265732 

25 0.365439093263182 0.0939902655253980 0.0813753254703597 0.0813753254703599 

35 0.341353356680001 0.0709936571240705 0.0582939643425133 0.0582939643425137 

50 0.323407733780370 0.0536207774061552 0.0408759738062723 0.0408759738062725 

Midpoint cubature for = cos (10a; -I- 20?/) ^—h cos(/?^ . 

Finally we consider the differentiable test function fz{x,y) = We obtain the following 

results 


Table 7 


N=M 


Error 

5 

1.790188901210552 

0.005 006 901 

10 

1.793908488488327 

0.001287 313 

20 

1.794870559302513 

0.000 325 242 

100 

1.795182719690259 

0.000 013 082 

200 

1.795192530240442 

0.000 003 271 


Midpoint cubature for 


True value is « 1.795195 802 05131 


If we compare the 2 x 9 x 35 = 630 valuations with precision 10“® in Table 3 with the present 
10000 evaluations with precision 10“^, we see how much better is our cubature formula, even if 
the integrand is a function. 
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4.3 Comparison with the generalized Peirce Rule 

S. De and K.J. Bathe discussed in [TH] (see also see [SU p. 65]) the following rule: For given N, let 
Pj, j = be the nodes of the Gauss quadrature Gn on with corresponding weights 

tPj,further a be a real parameter and M a natural number: then 




N M 


rPeirce,a / r\ r ( , - 2tT {s -\- a) - . 27r (s + a) 

Gr.. ’ (/) := [ y^cos- — -,VftSin 


3=1 


M 


is called the generalized Peirce rule. The rule of Peirce [41] is obtained by setting N 
M = 4m + 4, and a = 0. 


(40) 
m + 1, 


Remark 12 Let us note that the Discrete Polyharmonic Cubature for the constant weight function 
ru = 1 (hence K = 0 in h2(A) ) is identical with formula with a = 0. Thus the numerical 
evaluation of this formula can he carried out with our programme. 

The next table gives the results for the rule of Peirce for the first test function: 


Table 8 


N/M 

9 

25 

63 

83 

10 

6.49387212 

6.49387212 

6.49387212 

6.49387212 

15 

6.577936813 

6.577936813 

6.577936813 

6.577936813 

25 

6.647152541 

6.647152541 

6.647152541 

6.647152541 

35 

6.677370918 

6.677370918 

6.677370918 

6.677370918 

50 

6.700258414 

6.700258414 

6.700258414 

6.700258414 


Peirce Cubature rule for fow^^^ = (l + P + f —\- cos ip 


True value is —tt Ri 6.754 424 205 218 060 
20 

Even in the case N = 25, M = 25, with the number of evaluation points equal to NM = 25 • 25 = 
625, we have an error of 

6.754424205218060 - 6.647152541 = 0.107 271664 218 06 

In case of NM = 50 • 83 = 4150 evaluation points (which is about 6 times bigger than 625 ) the 
error is only about twice smaller: 

6.754424205218060 - 6.700258414 = 0.054165 791218 06 


The experiments with the test functions /i, /2 and /a have shown a similar behaviour. 

5 Experimental results for the weight function {x, y) 


In this section we discuss a weight function which is of quite different nature compared to the first 
weight function. The function 


(re*"^) := \y\ = |rsin(/j| 


(41) 


is simple in the sense that it is homogeneous in the variable r. However it has an infinite Fourier 
series since 


|sin(/:| = — 

TT 


4 

TT 


E 


fc=i 


cos {2kp) 
4fc2_l ■ 


The weight function is pseudo-definite since its orthonormalized Fourier coefficients have a 
definite sign: 


^"( 0 . 1 ) {-r) 



and u;( 2 fc.i) (r) 


4 1 

0F 4fc2- 1 


r for fc > 1. 
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We recall our main integration formula ((2) 


p 27 T pR oo ak pR 

Iw{f)= / / (re*‘^) w (re*"^) rdrd(p = / hkd) (^) (r) rdr. (42) 

Jo Jo k=oe=i 

5.1 Results for the discrete polyharmonic cubature 

Further we consider the test function 


/4 ix,y) = 30x^^. 


(43) 


Note that 

n 27r Q 

30r^^ cos^^ (ip) |sin y)| r'^drdip = — si 0.615 384 615 384 610. 

13 

Since {x, y) is a polynomial of degree 12 the Fourier coefficients satisfy fA,(k,t) (^) = 0 for fc > 13, 
hence, we may take K = 12. For M > 24 and > 12 the polynomial / {x,y) will be reproduced, 
i.e. holds 

<M.12) (/) = (/) = ■ 

This can be seen in the following table: 


Table 9 


N\M 

9 

25 

63 

83 

10 

0.5609353695139790 

0.6153846153846160 

0.6153846153846160 

0.6153846153846160 

15 

0.5609353656165750 

0.6153846153846150 

0.6153846153846150 

0.6153846153846150 

25 

0.5609353655541600 

0.6153846153846170 

0.6153846153846170 

0.6153846153846170 

35 

0.5609353655539850 

0.6153846153846130 

0.6153846153846140 

0.6153846153846140 

50 

0.5609353655539860 

0.6153846153846170 

0.6153846153846170 

0.6153846153846170 

Error 





N\M 

9 

25 

63 

83 

10 

0.0544492458706369 

0.0000000000000000 

0.0000000000000000 

0.0000000000000000 

15 

0.0544492497680410 

0.0000000000000010 

0.0000000000000010 

0.0000000000000010 

25 

0.0544492498304560 

0.0000000000000010 

0.0000000000000010 

0.0000000000000010 

35 

0.0544492498306309 

0.0000000000000030 

0.0000000000000020 

0.0000000000000020 

50 

0.0544492498306299 

0.0000000000000010 

0.0000000000000010 

0.0000000000000010 


Discrete Polyharmonic cubature for K = 12 and = 30x^^ \y\ 

True value is ft: 0.615 384 615 384 616 


Next we consider the test function 


f 5 {x,y) :=\y\. (44) 

Since it is not smooth we should not expect to obtain too good results for the discrete polyharmonic 
cubature. The exact value of the integral can be computed: 

C 2 C C’' 2 

h (|y| |y|) = / / |rsin(p| diprdr = / r'^dr- / |sin(p| d(p 

Jo Jo Jo Jo 

= irr ft 0.785 398 163 397448 
4 

We took K = 12, as in the last experiment, and we obtained a table of results where the error is 
almost the same for all N and M, and is around 10“^. If we take K = 22 the results are much 



















better: 


Table 10 


N\M 

9 

25 

63 

83 

10 

0.785206660 

0.785352337 

0.785367124 

0.785369362 

15 

0.785208297 

0.785358970 

0.785373081 

0.785375274 

25 

0.785208235 

0.785361119 

0.785374994 

0.785377171 

35 

0.785208149 

0.785361440 

0.785375276 

0.785377452 

50 

0.785208109 

0.785361541 

0.785375364 

0.785377539 

Error 





N\M 

9 

25 

63 

83 

10 

0.0001915037 

0.0000458267 

0.0000310393 

0.0000288009 

15 

0.0001898666 

0.0000391939 

0.0000250824 

0.0000228890 

25 

0.0001899281 

0.0000370443 

0.0000231697 

0.0000209919 

35 

0.0001900142 

0.0000367231 

0.0000228870 

0.0000207116 

50 

0.0001900544 

0.0000366227 

0.0000227993 

0.0000206248 


Discrete Polyharmonic cubature for K = 22 and = \yf 


True value is « 0.785 398163 397448 


Finally, we look again at the oscillating test function ([36|) . namely, /2 {x, y) = cos (lOx + 20?/) .Using 
the expansion 


('-lU 

C°-'^ (x, y) := cos (ax + by) = ^ 72^ 
one can prove that the Fourier coefficients C^ 2 k+i i) (^) 

^( 2 M) = 27ry2fc.^ {‘Pa,b) (-1)^ J 2 k (r\/ (a^ + 6^)^ 


where J 2 k are the Bessel functions, see formula (l38l) . and (pa,b is the angle given by ta,Tnpa,b = 
It follows that 


(/) = 4 J Jo (rV a? + 62^ r'^dr + ^ 4 fe 2 °!.^x (r\/a2"+^) r^dr. (45) 

We may estimate the error using this expression. 

In formula m, for the Discrete Polyharmonic Cubature, we have taken K = 22. We 
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obtain the following results: 


Table 11 


N\M 

9 

25 

63 

83 

10 

0.096718846427391 

0.014472433304185 

0.014477271351135 

0.014477271351135 

15 

0.096642162991824 

0.014472441635349 

0.014477279682299 

0.014477279682299 

25 

0.096593708566055 

0.014472441635349 

0.014477279682299 

0.014477279682299 

35 

0.096592741482745 

0.014472441635350 

0.014477279682299 

0.014477279682299 

50 

0.096597855764522 

0.014472441635349 

0.014477279682299 

0.014477279682299 

Error 





N\M 

9 

25 

63 

83 

10 

0.082241566745092 

0.000004846378115 

0.000000008331165 

0.000000008331165 

15 

0.082164883309524 

0.000004838046951 

0.000000000000001 

0.000000000000001 

25 

0.082116428883756 

0.000004838046950 

0.000000000000000 

0.000000000000000 

35 

0.082115461800446 

0.000004838046950 

0.000000000000000 

0.000000000000000 

50 

0.082120576082222 

0.000004838046950 

0.000000000000000 

0.000000000000000 


Discrete Polyharmonic cubature for K = 22 and = cos(10rcos:/? + 20rsin(/3) |rsin(p| 


True value is 0.0144772796822995 

5.2 Comparison with piece-wise midpoint rule 

We want to give briefly the results for the piece-wise midpoint rule for the test function (cc, y) = 
30a;^^ with the weight {x,y) = \y\. 


rmid 

^ 5,5 

(30x12 

I 2 /I) 

= 0.173 

359 

053 

300 

102 

rmid 
'^10,10 

(30x12 

I 2 /I) 

= 0.795 

909 

972 

453 

979 

rmid 

-^20,20 

(30x12 

\y\) 

= 0.640 

355 

591 

783 

074 

rmid 

-^ 40,40 

(30x12 

\y\) 

= 0.620 

920 

690 

132 

442 

rmid 

-^100,100 

(30x12 

I 2 /I) 

= 0.616 

243 

873 

682 

115 

rmid 

-^200,200 

(30x12 

I 2 /I) 

= 0.615 

598 

519 

244 

782 

rmid 

-^ 500,500 

(30x12 

I 2 /I) 

= 0.615 

418 

799 

448 

66 

True value is 

« 0.615 

384 

615 

384 

616 


We see that even with 500 x 500 evaluation points the error is around 0.0001. On the other hand, 
for the test function /s {x,y) = \y\ the piecewise midpoint rule is very good, since the integrand is 
just \y\ but we do not have exactness. 

We omit a discussion of the Peirce rules since the results are very similar to the case of the 
midpoint cubature rule. 

6 Aspects of the Numerical Implementation 

The discrete polyharmonic cubature (l20)l is defined for pseudo-definite weight functions w (re*"^) 
on the disk. It is important to determine numerically the nodes and weights for the quadrature of 
degree N for the univariate integrals 

^ L {y/p) dp 

with high accuracy. The nodes are the zeros of the orthogonal polynomial (p) of degree N with 
respect to the measure {^/p) ■ Thus any implementation of our algorithm depends on 

reliable and stable software for finding the nodes and weights of the corresponding quadrature. We 





















have judiciously chosen the weights and such that 

W(k,e) iVp) = (1 - 

are weight functions of Jacobi type, i.e. of the form a;“ (1 — x)^ with a, /3 > — 1 in the interval [0,1]. 
For these weight functions fast and highly accurate programmes for finding the nodes and weights 
are available (see e.g. [22] or the website: https://www.cs.purdue.edu/archives/2002/wxg/codes) 
which we have used in our Matlab programs. A more general situation, namely when the Fourier 
coefficients W(k^t) (y^) linear combinations of Jacobi type weight functions, is easy to handle. 
In the case that w^k^n) {^/p) is a general pseudo-definite function the user has to search or develop 
numerical procedures for finding the nodes and coefficients with high accuracy. Another idea, which 
we did not pursue in detail, is to approximate W(^k,e) (\/p) by means of a Bernstein polynomial 
~ pY of some degree d, which is then a linear combination of Jacobi weights. 
However, a problem might be the approximation rate for d —> oo, which is notoriously not very 
good. 

In our Matlab program we used the readily implemented Matlab function for the Fast Fourier 
Transform: then the Fourier series is in the form 

OO 

w (re*'^) = ^ Wk (r) 

k— — oo 


with the usual relations wq (r) = ^W(^o,i) (f) and 

W(k,i) {r) = ^ (wk (r) + Wk (r)^ and wi^k, 2 ) {r) = ^ (wk (r) - Wk (r)) . 

Here Wk (r) := ^ w (re®'^) For the computation of the approximation to the 

Fourier coefficients f(k,e) (j) of the integration function / (a:, y) we use the Fourier series represen¬ 
tation 

oo .. „ 27 r 

f{re^^)= ^ /, where/fe(r):=—y f [re^^] d^. 


Hence, the integral (HU can be written in the form 

oo pR 


oo .R oo .R 

i(f)= E / 27r/fc (r) w-k (r) rdr = '^ / nfk {y/p) W-k {y/p) dp. 

fe=-oo'^0 fe=-oo‘^0 


(46) 


In our program we choose always odd M, and the discrete Fourier approximation to fk (r) given 

by 


M 


(r) := (r) - (^) = e 


— 2772^ 


S = 1 


which is the link to the Fast Fourier transform. By subdividing the interval [0, 27r] into M subin¬ 
tervals we obtain an approximation 






^—‘I'KiksjM. 


of fk (r) which is just the Discrete Fourier transform (DFT) for the data points / for 

s = I,...,M. 




7 Concluding Remarks 

The discrete polyharmonic cubature formula (/) provides excellent numerical results for 

integrating functions on the disk in the plane with respect to a weight function w. A possible 
drawback for applications might be the high number of evaluations points needed in the formula 
given by {2K — 1) ■ N ■ M evaluation points. Since the evaluation of a function value / (x) might 
be very costly, it is a natural question whether one could modify the formula so that we would 
need less function evaluations. 

In a forthcoming paper [32] we introduce a new cubature formula related to the above formula 
(1211) with a high degree of ” approximative exactness”, which uses a spline approximation in direc¬ 
tion r of the function (x). The point evaluations for this formula are on a regular grid in both 
directions tf and r, and the number of knots is Ni x M, where Ni is the number of spline knots 
in direction r. We call this a hybrid cubature since we combine spline methods with a cubature 
formula. In mathematical terms the hybrid formula is of the form: 


, if afc AT 

fe=o e=i0=1 

{t) is a univariate spline interpolation function with nodes 
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where SPL 


{/s; (^™)} 


Ni 

m—O 

N, 


for the data |/(^) 


{/(M) X (^7) 
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